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ABSTRACT 

Resistivity and viscosity have a significant role in estabhshing the energy 
levels in turbulence driven by the magnetorotational instability (MRI) in local 
astrophysical disk models. This study uses the Athena code to characterize the 
effects of a constant shear viscosity u and Ohmic resistivity r] in unstratified 
shearing box simulations with a net toroidal magnetic fiux. A previous study of 
shearing boxes with zero net magnetic field performed with the ZEUS code found 
that turbulence dies out for values of the magnetic Prandtl number, Pm = ^ hi 
below Pin ~ 1; for Pm > 1, time- and volume- averaged stress levels increase 
with Pm. We repeat these experiments with Athena and obtain consistent re- 
sults. Next, the infiuence of viscosity and resistivity on the toroidal field MRI 
is investigated both for linear growth and for fully-developed turbulence. In the 
linear regime, a sufficiently large v oy r] can prevent MRI growth; Pm itself has 
little direct infiuence on growth from linear perturbations. By applying a range 
of values for v and rj to an initial state consisting of fully developed turbulence 
in the presence of a background toroidal field, we investigate their effects in the 
fully nonlinear system. Here, increased viscosity enhances the turbulence, and 
the turbulence decays only if the resistivity is above a critical value; turbulence 
can be sustained even when P^ < 1, in contrast to the zero net field model. 
While we find preliminary evidence that the stress converges to a small range of 
values when v and rj become small enough, the infiuence of dissipation terms on 
MRI-driven turbulence for relatively large rj and v is significant, independent of 
field geometry. 

Subject headings: accretion, accretion disks - black hole physics - (magnetohy- 
drodynamics:) MHD 
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Introduction 



Disk accretion is a fundamental process of many astrophysical phenomena, from nearby 
young stellar objects to immensely luminous distant quasars. Understanding the mechanism 
for removing angular momentum from a fluid element, thereby allowing accretion to occur, 
is essential to understanding these s ystems. Orbiting, magnetize d gas is unstable to the 
magnetorotational instability (MRI) (IBalbus fc Hawleylll99ll . Il998l ): all that is required is a 
subthermal magnetic field sufficiently coupled to differentially rotating gas with a negative 
outward angular velocity gradient. The MRI leads to turbulent flow, resulting in Maxwell 
and Reynolds stresses that efficiently transport angular momentum and drive accretion. 
However, it is still uncertain what determines the amplitude of the magnetic energy and 
stress in saturated MRI-driven turbulence. 

Because linear analysis can offer only limited guidance, numerical simulations have been 
used to investigate the properties of MRI-driven turbulence. Most simulations of the MRI 
employ the shearing box approximation, in which the simulation domain consists of a local 
corotating patch of accretion disk, small enoug h to expand the MH D equations into Cartesian 
coordinates and ignore curvature terms (see iHawley et al.lll995l ). This approximation, in 
its simplest form, reduces the problem to its basic ingredients: differential rotation and 
magnetized fluid. It is hoped that a more complete understanding of this simple model 
will provide insights into the mechanisms that determine the stress levels in astrophysical 
systems. 



The flrst shearing box simulations (IHawley et al.l Il995l . Il996l ) found that the presence 
of a net magnetic fleld and its orientation play a role in setting the amplitude of the MRI 
turbulence. A net magnetic fleld results from currents located outside of the computational 
domain and cannot change as a result of the evolution. For zero net fleld simulations, on 
the other hand, complete decay of the fleld is possible. Net vertical flelds gave the largest 
turbulent energies, with the energy level approximately proportional to the background fleld. 
Net toroidal flelds behave similarly, but with a smaller energy for the same background fleld 
strength. Zero net fleld simulations saturated at levels comparable to those seen in net 
toroidal fleld cases. 

In subsequent years, there have been many shearing box simulations which conflrm these 
qualitative behaviors, but the factors that determine the amplitude of the turbulent energy 
still remain uncertain. Some studies have examined numerical effects such as computational 
domain size and resolution, and others have looked at physical parameters such as back- 
ground fleld st rength and gas pr essure. An investigation of the influence of gas pressure 
carried out by ISano et al.l (120041 ). for example, found an extremely weak pressure depen- 
dence. Even here, the influence of the gas pressure depends on the magnetic fleld geometry. 
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Blackman et al.l (120081 ) examined the results of an ensemble of shearing box simulations 
taken from the literature and found that a(3 is generally constant, where a is the total stress 
divided by the gas pressure, and (3 is the ratio of thermal to magnetic pressure. In other 
words, stress is proportional to magnetic energy. 

A physical influence that has, until recently, received less attention is physical dissipa- 
tion, namely shear viscosity v and Ohmic resistivity rj. The linear dispersion re l ation for the 
vertical field MRI in the presence of z/ and r] was derived by iBalbus fc Hawleyl (Il998l ) . Both 
terms can reduce the effectiveness of the MRI. In the linear regime, viscosity damps the MRI 
growth rates and changes the wavelength of the fastest growing mode, but does not alter the 
wavenumbers that are unstable. Resistivity introduces a cutoff on the unstable wavelengths 
where the res istive diffusion time be comes comparable to the MRI growth time (see, e.g., the 
discussion in Masada fc Sanol 20081). Nonaxisyrn metric MRI modes with Ohmic resistivity 
were examined by iPapaloizou fc TerquemI (119971 ) . They found that resistivity reduces the 
amplification of such modes, and if large enough, can stabilize the toroidal field MRI. 



Si mulations bvlHawlev et a. 



mm . IZieder k Riidigeil (120011 ) . and 



(ll996).ISano et al.l (1998h [Fleming et al.l (l2000h . lSano k Inutsuka 



Sano k Stond (120021 ) have investigated the influence 



of a nonzero Ohmic resistivity on the saturated state. The main result of these studies is 
that increasing the resistivity leads to a decrease in turbulence, independent of the magnetic 
field configuration. In zero net field models, the effect of r esistivity on the tur bulence is 
larger than one might expect from the linear MRI relation (jFleming et al.ll2000l ). with the 
turbulence decaying to zero for relatively low values of resistivity. 



Recently, the work of iFromang et al.l (120071 ) (hereafter F07) and iLesur k Longaretti 



(120071 ) has sparked new interest in the effects of non-ideal MHD on the MRI. F07 showed 
that both resistivity and viscosity are impo rtant in determining the s tress level in MRI 
turbulent flows with zero net magnetic field. iLesur k Longarettil (120071 ) came to the same 
conclusion for models with a net vertical field. The results were characterized in terms of the 
magnetic Prandtl number, defined as Pm = ^/v- In these simulations, the saturation level 
increases with increasing P^. F07 also find that for the zero net field case, there exists a 
Pjn below which the turbulence dies out, and that this critical Pm decreases with decreasing 
viscosity (at least for the range in viscosity and resistivity examined in the paper). 

One magnetic field geometry that has not yet been explored with both physical resis- 
tivity and viscosity is that of a net toroidal field . Such fields could be the most relevant to 
astrophysical disks. Following the arguments of iGuan et al.l (120091 ) and references therein, 
both global and local disk simulations as well as observations of disk galaxies show a dom- 
inance of toroidal field over other field components. Indeed, the background shear flow 
naturally creates toroidal field from radial field. It seems likely that any given region of an 
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accretion disk will contain some net azimuthal field. 

In this paper, we perform the first investigation of the toroidal field MRI in the presence 
of both viscosity and resistivity and compare the outcomes with those obtained for zero net 
and net vertical field simulations. The structure of the paper is as follows. In § |2l we describe 
our algorithm, parameters, and tests of our viscosity and resistivity implementation. For 
comparison purposes, we reexamine the simulations of F07 with our code in § [31 Our main 
results, focusing on the toroidal field simulations, are presented in §|H We wrap up with our 
discussion and conclusions in § O 



2. Numerical Simulations 



In this study, we use the Athena code, a second-order accurate Godunov scheme for 
solving the equations of MHD in con servativ e forrn using the dimensionally unsplit corner 
transport upwind (CTU) method of IColella couple d with the third-order in space 

piecewise parabolic method (PPM ) of IColella fc WoodwardI (jl98J) and a constrained trans- 
port (CT; lEvans &: Hawleylll988l ) algorithm for preserving the V • .B = constra int. A de- 



scription of this al g orithm and various test proble ms is given in [Gardiner &: Stond (l2005bl ) , 
Gardiner &: Stond (120081 ) , and IStone et al.l (120081 ) . For the present study, we have added 
physical dissipation in the form of a constant kinematic shear viscosity and Ohmic resistiv- 
ity using operator splitting, as described in more detail below. Bulk viscosity is ignored. 

The shearing box approximation is a model for a local region of a disk orbiting at a 
radius R whose size is small compared to this radiu s, allowing us t o exp and the equations of 
motion in Cartesian form, as described in detail by lHawley et al.l (119951 ). The box corotates 
with an angular velocity Q corresponding to the value at the center of the box. The shearing 
box evolution equations with viscosity and resistivity are given by iBalbus &: Hawleyl (119981 ) 
and are: 

rin 

(pv) = 



1 



dpv 
~dt 



1 



+ V-{pvv-BB) + V P + -B^ = 2qpn^x -2n X pv + V ■{piyVv)+V -puW ■ v 



dB 

'dt 



V X {v X B - r]V X B) 



(2) 
(3) 



where p is the mass density, pv is the momentum density, B is the magnetic field, P is the 
gas pressure, and q is the shear parameter, defined as g = — dlnfi/dlni?. We use q = 3/2, 
appropriate for a Keplerian disk. We assume an isothermal equation of state P = pc^, where 
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Cs is the isothermal sound speed. Shear viscosity and Ohmic resistivity are denoted by v 
and T] respectively. Note that our system of units has the magnetic permeability = 1. The 
first source term on the right-hand side of equation ([2]) corresponds to tidal forces (gravity 
and centrifugal) in the corotating frame. The second source term in equation ([2]) is the 
Coriolis force. Finally, we have omitted the vertical component of gravity, making these 
"unstratified" shearing box simulations. 

Adapting the Athena code to the shearing box problem requires adding the tidal and 
Coriolis force terms and implementing the shearing-periodic boundary conditions at the 
X boundaries. The source terms are included in the algorithm in a directionally unsplit 
m anner, consistent with the CTU algorithm. We do not use the Crank-Nicholson method 



of iGardiner &: Stond (l2005a( ) that ensures precise conservation of epicyclic energy. We have 



found this a dded complexity to be unnecessary for simulations dominated by the MRI (see 
arguments in Simon et al.ll2009l) . The shearing-periodic boundary conditions are described in 
Simon et al.l (120091 ) . Quantities are linearly reconstructed in the ghost zones from appropriate 
zones in the physical domain that have been shifted along y to account for the shear across 
the boundaries. Furthermore, the y momentum is adjusted to account for the shear across 
the X boundaries as fluid moves out one boundary and enters at the other. 

Note that to preserve a quantity to machine precision across a grid boundary such as 
the shearing-periodic boundary (or a boundary between different grids in a mesh-refinement 
scheme), it is necessary to reconstruct a quantity's flux (or for the magnetic fi eld, the elec- 



trom otive force, EMF) at the boundary, rather than the quantity itself (see ISimon et al. 



20091 ) . To conserve total vertical field flux, for example, we reconstruct the y EMF at the 
X boundaries. This is essential, given the strong effect that a net vertical field has on the 
turbulence level. The perfect conservation of net toroidal flux is not as important, and as 
ensuring its precise conservation involves a more complex procedure, we allow the By flux 
to be conserved only to truncation level. Note that in our simulations initialized with a net 
toroidal fleld, this truncation error results in a loss of net By flux from the domain; ~ 5-10 



of the initial toroidal fleld is lost per 100 orbits for our high resolution, sustained turbulence 
simulations. This corresponds to a background value of ~ 110-120 at 100 orbits. While 
this truncation error does not appear to have any signiflcant affect on the turbulent energy 
levels in our simulations, it may become important to conserve By to roundoff level for longer 
evolution times. The radial flux, B^-, will automatically be conserved to machine precision 
because its evolution is determined by EMFs on the periodic y and z boundaries. We also 
reconstruct the density flux on the x boundaries to conserve the total mass in the domain to 
machine precision. The systematic difference between the calculation of outwar d and inward 



fluxe s at the shearing x boundaries can lead to a loss of mass from the grid (ISimon et al. 



20091 ). We do not reconstruct momentum fluxes at the boundaries as the source terms will 
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prevent roundoff level conservation of momentum. 

Both the viscosity and resistivity are added via operator splitting; the fluid variables 
updated from the CTU integrator are used to calculate the viscous and resistive terms on the 
right-hand side of equations ([2]) and ([3]). These terms are discretized in a flux-conservative 
manner consistent with the Athena algorithm. In particular, the third and fourth terms on 
the right-hand side of equation ([2]) are written so that puVv and (l/3)pz/V ■ v are defined 
as fluxes at the cell faces. Taking the divergence of the third term and the gradient of the 
fourth term via finite-differencing ensures that momentum conservation is not violated by 
the viscous terms. The resistive contribution to the induction equation is added in a manner 
consistent with the EMFs; the term 77V x S is computed at cell corners to ensure that 
when differenced via the curl operator, V ■ S = is maintained. Note that this resistive 
contribution to the EMF must also be reconstructed at the shearing-periodic boundaries in 
order to preserve Bz precisely. 

More generally, the viscous term in equation ([2]) can be written in a flux-conservative 
manner as V • T where T is a viscous stress tensor defined as 

^ / dvi dvj 2 ^ ^ \ , ^ , 



where the indices refer to the spatial components (ILandau fc Lifshitzl Il959l ). For simplic- 
ity, we have used the form as in equation ([2]) which is equivalent to equation (jl]) assuming 
that pv is spatially constant. We have performed a few shearing box experiments with 
both implementations, and find no significant differences in turbulent stress evolution. In 
particular, we restarted a few simulations using the form in equation (jlj). We found that 
the volume-averaged magnetic energies are initially indistinguishable between the two ap- 
proaches. Due to the chaotic nature of the MRI, the two curves eventually diverge, but 
nevertheless maintain the same time average. 

The addition of viscosity and resistivity places an additional constraint on the time step. 
At = C,MIN (^AtcTu, 0-75^, ^-^^^ ' (5) 

where Co is the CFL num ber (Cn = 0.4 here ), AtcTU is the time step limit from the main 
integration algorithm (see Istone et allboosh . and A is the minimum grid spacing, A = 



MIN{Ax, Ay, Az). Several three-dimensional tests of viscosity and resistivity revealed that 
if the viscous or resistive time step is close to AtcTU, the evolution becomes numerically 
unstable. This problem was remedied by multiplying the viscous and resistive time steps 
by 0.75. The additional 4/3 factor in the denominator of the viscous time step results from 
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the last term on the right-hand side of equation ([2]). This can be most easily understood 
by considering a one-dimensional problem, in which case the effective v value increases by a 
factor of 4/3 due to the compressibility term. Therefore, the effective v that goes into the 
time step calculation is taken as (4/3)z/. Note that most of our simulations will have v and 
r] sufficiently small that the viscous and resistive time steps are large compared to AtcTU- 
In fact, only the simulations with the largest values of ri and v reach the diffusion limit on 
At. 



2.1. Tests of Physical Dissipation 



We performed a number of problems to test the implementation of viscosity and resistiv- 
ity within Athena. Resistivity was tested by solving the diffusion of a current sheet along one 
dimension; a uniform magnetic field is initialized with a change in sign acr oss one grid zone. 
This problem has a simple analytic solution (see e.g., iKomissarovl 120071 ). The agreement 
between the numerical and analytic solution was excellent. By replacing the magnetic field 
with a uniform velocity flow, the identical test can be performed for the viscosity. Again, 
the numerical solution agreed with the analytic solution. 

Next, we initialized a uniform vertical magnetic field in a shearing box with nonzero 
viscosity and resistivity and measured the growth of various MRI modes in the linear 
regime. We cornpared the measured values with those from analytic linear theory (see e.g., 
Masada fc Sandl2008l : iPessah fc Charul2008l ) and found excellent agreement for a wide range 
in viscosity and resistivity. 

Finally, we examined the propagation of small amplitude, isothermal sound and Alfven 
waves in the presence of viscosity and resistivity. Again, the numerical solution can be 
compared directly to an analytic solution. These tests were done in one, two, and three 
dimensions; in the multidimensional tests, the propagation direction of the wave was chosen 
to be along the grid diagonal. The resistivity was tested via the decay of the Alfven waves, 
and the viscosity was tested via the decay of the sound waves. The error as a function of x 
resolution for two of these tests is given in Fig. [H The error is calculated from the square 
root of the sum of the squared errors in the density and momenta (for the sound wave) and 
the density, momenta, and magnetic field (for the Alfven wave). The solution to each wave 
converges at a rate very close to second order, shown by the dashed line. 
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2.2. Sheeiring Box Pcirameters 

The shearing box used in this study has radial size = 1, azimuthal size Ly = 4, and 
vertical size Lz = 1. Most of the simulations presented here use 128 x 200 x 128 equally 
spaced grid zones; some simulations use half the number of zones in each direction. The 
initial velocity is v — —qVlxy, with q = 3/2, Vt — 0.001, and —L^jl < x < The 
isothermal sound speed is Cg = 0,H where H is the scale height. With — H , we have 
Cg = Lzi^, and with p—1, the initial pressure is P = P^^L^ — 10~^. 

The dissipation terms v and rj are parameterized in terms of the Reynolds number. 



Re = (6) 



the magnetic Reynolds number. 



Rm = (7) 



and the magnetic Prandtl number. 



P.^-^ = ^. (8) 

Tj Re 

Since the properties of the MRI are more directly determined by the Alfven speed rather 
than the sound speed, another useful dimensionless quantity is the Elsasser number. 



where va is the Alfven speed. With Cg = and P — 2cl/v\, we can relate Rm to A, 



A = -Rm. (10) 

In addition to the exphcit dissipation terms, there will also be some effective diffusion 
due to numerical effects. Generally speaking, numerical diffusion will not behave in the 
same manner as physical diffusion (e.g., it is not a simple function of a gradient in field or 
velocity); numerical diffusion generally has a much stronger effect at small scales than at large 
scales. Also the effects of numerical diffusion may be different from one type of simulation to 
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another. By calculating numerical losses at high wavenm bers in Fouri e r spac e and modeling 
those as if they were physical viscosity and resistivity, ISimon et al.l (120091 ) quantified the 
numerical dissipation of Athena. They found that the effective Rm for the zero net field and 
net z field simulations at N;^. = 128 were 20000 and 8000 respectively, and 7000 and 5000 
for Nrc = 64. The effective Pm is ~ 2 for these simulations. Since numerical dissipation is 
problem-dependent, these numbers should be regarded as estimates, and their values will 
likely be somewhat different in different applications. However, they serve as a guideline 
for including physical dissipation. In the present study, numerical and physical dissipation 
may be comparable at large wavenumbers for Re, Rm > 10000. The physical dissipation in 
some of our simulations may fall into this marginally resolved regime. Nevertheless, we can 
explore a large enough range in Re and Rm values to observe clear effects due to viscosity 
and resistivity. 



Zero Net Flux Simulations 



Fromang fc Papaloizoij (120071 ) and iPessah et al.l (120071 ) presented the surprising result 
that for zero net field shearing box simulations without any explicit di ssipation terms, the 
steady-state turbulent energy decreases with increasing grid resolution. ISimon et al.l (120091 ) 
obtained the same result for zero net field simulations without explicit dissipation using the 
Athena code. These results pointed to the importance of including explicit dissipation terms 
in such simulations. 

F07 showed that turbulent activity is strongly infiuenced by these dissipation terms; the 
saturated stress increases with increasing P^. Here we return to the zero net field problem 
and include the dissipative terms to compare with the results of F07. The simulations are 
initialized with B = ^/2P/psm[{2^T/L^)x]z where (3 = 400. These runs are labeled SZ for 
sinusoidal z-field and have resolution A^^ = 128, Ny = 200, A^^ = 128. The viscosity and 
resistivity in these simulations are chosen to reproduce the calculations of F07. The initial 
state is perturbed in each grid zone with random fiuctuations in p at amplitude Sp/p = 0.01. 
The SZ simulations are listed in Table [H The column labeled "Turbulence?" states whether 
or not the turbulence was sustained in a given simulation. The column labeled "a" gives 
the resulting turbulent stress in terms of the dimensionless value a = {{pv^Svy — B^By)) / Po, 
with 6vy = Vy + qflx. Po is the initial gas pressure and the double bracket denotes a time and 
volume average. The time average is calculated from orbit 20 until the end of the simulation, 
and as is the case throughout this paper, volume average refers to an average over the entire 
simulation domain. 



The results of these simulations are consistent with those of F07. For example, F07 lists 
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a values for a Re = 3125 and Rm = 12500 model run with four different codes, including 
ZEUS. These values range from a = 0.0091 to 0.011; we obtain 0.013. The increase in 
turbulent energy levels with is demonstrated by a series of simulations with the same Rm 
and increasing viscosity. For example, for a constant Rm ^ 12800 (some of the simulations 
had Rm = 12800 while others had Rm = 12500; see F07), Pm values were varied by factors 
of 2 from 1 to 16. Sustained turbulence was seen for Pm > 4 with a values increasing from 
0.0091 for Pjn = 4 to 0.019 and 0.044 for Pm = 8 and 16 respectively. The Athena runs have 
a values of 0.013, 0.026, and 0.046. These data are plotted in Fig. [21 which shows that the 
increase in a with Pm is nearly linear. 

The largest differences between the Athena simulations and the ZEUS simulations of F07 
lie in the marginally turbulent cases. For example, we find decaying turbulence for Re = 1600, 
Pm = 4, whereas ZEUS produces sustained turbulence for these parameters. Figure [3] shows 
the volume-averaged magnetic energy density normalized by the gas pressure versus time for 
the three Pm values at Re = 1600. The lowest Pm simulation decays quite rapidly, whereas 
the Pm = 4 case takes roughly 60 orbits to decay. Differences in the numerical properties of 
Athena and ZEUS might account for these results, given the sensitivity to numerical factors 
as shown by zero net field simulations. We also note that we use a slightly larger domain size 
in y than in F07. The boundary in parameter space between sustained turbulence and decay 
is unlikely to be hard and fast, and detailed numerical surveys that attempt to define that 
boundary are probably not warranted. Some such studies may, however, provide additional 
insights into the sensitivity of the MRI turbulence to specific numerical factors. 



Table 1. Zero Net Flux Simulations 



Label 


Re 


Pm 


Rm 


Turbulence? 


a 


SZRe800Pm4 


800 


4 


3200 


No 




SZReSOOPmS 


800 


8 


6400 


Yes 


0.031 


SZRe800Pml6 


800 


16 


12800 


Yes 


0.046 


SZRel600Pm2 


1600 


2 


3200 


No 




SZRel600Pm4 


1600 


4 


6400 


No 




SZRel600Pm8 


1600 


8 


12800 


Yes 


0.026 


SZRe3125Pml 


3125 


1 


3125 


No 




SZRe3125Pm2 


3125 


2 


6250 


No 




SZRe3125Pm4 


3125 


4 


12500 


Yes 


0.013 
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4. Toroidal Field Simulations 

To examine the effect of viscosity and resistivity on the MRI with a net toroidal field, 
we have run a series of simulations initialized with B = a/2P/ (3y, where (3 = 100, and 
with varied Re and Rm values. Re ranges from 100 to 25600, and ranges from 0.25 to 
16 (though, in some simulations, we set either rj or u equal to zero). We will consider the 
influence of the physical dissipation terms on two types of problems: the linear MRI growth 
regime, and fully nonlinear turbulence. 



4.1. The Linear Regime 



The linear nonaxisymmetric MRI was first examined by lBalbus &: Hawleyl (Il992l ). For 
nonaxisymmetric modes, the MRI tends to be most robust in the presence of a poloidal 
fiel d. However, even the p urely toroidal field case is unstable, athough, as emphasized 
by iBalbus fc HawlevI ([1992), that case is somewhat singular. As always with the ideal 
MRI, the most unstable mode has k ■ va — ^- The linear analysis is complicated by the 
background shear which causes radial wavenumbers to evolve with time. Amplification of a 
given mode occurs when the wavenumber ratio k/k^ goes through a minimum as the radial 
wavenumber swings from leading to trailing. In general, the purely toroidal MRI favors high 
kz wavenumbers and small values of ky/kz, in contrast to the vertical field MRI where the 
wavenumber kz of the most unstable mode is determined by the Alfven speed. 



Papaloizou fc TerquemI (119971 ) examined the toroidal field MRI with the addition of 



resistivity. They point out that because kx grows arbitrarily large, all linear modes will 
eventually damp out in the presence of resistivity. For small enough resistivities, however, 
there can be a period of growth when A;^ ~ 0. For the MRI to become self-sustaining, this 
growth has to continue long enough for the perturbations to reach nonlinear amplitudes. 
Resistivity is also particularly important for th e pure toroidal field MRI bec ause large kz 
is favored for mode growth. Equation (32) of iPapaloizou &: TerquemI (119971 ) provides an 
approximate condition for transient amplification of the MRI in the presence of resistivity. 
For Keplerian shear and for modes where k ■ va ~ ^, this reduces to the condition 



k^ri ~ n. 



(11) 



In other words, there is no amplification of modes for which the diffusion time is comparable 
to the orbital frequency. Although viscosity was not included in the analysis, one might 
expect it to be similarly influential. 



Simulations of the linear growth of the MRI in the presence of resistivity for a purely 
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toroidal (3 = 100 initial field were first carried out by lFleming et al.l (|2000| ) using a ZEUS code 
with an adiabatic equation of state. For this field strength, the critical MRI wavelength in 
the azimuthal direction is 2itva/Q ~ H. They found field decay for a Rm = 2000 simulation, 
but field growth to turbulent saturation for Rm = 5000 and above. 

In this section, we follow the growth of the MRI in a shearing box with a purely toroidal 
field while including both resistivity and viscosity. The system is seeded within each grid 
zone with random perturbations in p at amplitude Sp/p = 0.01. The simulations were run 
at two resolutions, A^^^ = 64, Ny = 100, = 64 and = 128, A^'j^ = 200, = 128 and are 
labelled YL for y-field, linear regime. In this standard set of simulations, the range of Re 
examined runs from 800 to 25600, and the range of Rm is from 400 to 102400. Table [2] lists 
these simulations. The last two columns state whether or not MRI growth is observed for 
the N;j. = 64 and A^^ = 128 resolutions, respectively. A dash in either of these columns means 
that the simulation was not run at that particular resolution. MRI growth is defined by the 
evolution of the volume-averaged magnetic and kinetic energy components. A simulation 
is considered to have zero growth if after 20-40 orbits, the various energy components are 
either decaying or constant in time without any indication of exponential increase. Growth 
to saturation is observed in cases when Re and Rm are at 6400 and above. 

Clearly, a sufficiently large viscosity or resistivity can inhibit growth. But what about 
the very high or very low Pm limits? To approach that question, we carried out simulations 
where only u oi t] was nonzero. These experiments were done at the A^^,' = 64 resolution. In 
our first experiments, we set t] to zero and Re to 100 and 800. The Re = 800 run showed 
growth to saturation, but the Re = 100 case had no growth. Next we set u to zero and Rm 
to 800 and 1600. The lower resistivity {Rm = 1600) grew to saturation, whereas the higher 
resistivity {Rm = 8 00) did not. A lt hough the existence of a critical Rm value is consistent 



with the results of [Fleming et al.l (120001 ). the value of Rm at which growth is prevented 
is smaller here than what they found. We note that there remains unavoidable numerical 
dissipation associated with grid scale effects, which will make the value of a critical Rm 
obtained through simulations somewhat dependent on algorithm and resolution. 

The effect of numerical resolution is not necessarily obvious. Consider model YLRe3200Pm2, 
which has Re = 3200 and Rm = 6400, and model YLRe6400Pm0.5, which has these values 
reversed. In both cases, the A^^ = 64 simulations show growth but the A^^ = 128 models 
do not. One difference between the two resolutions is in the initial perturbations. While 
the density perturbations have the same amplitude in both resolutions, the higher resolution 
initial density is given power at smaller scales because the perturbations are apphed to each 
grid zone. This leads to a smaller amplitude for each Fourier mode. Does this account for the 
difference seen in these two resolutions? To investigate this, we ran both A^^; = 64 versions of 
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YLRe3200Pm2 and YLRe6400Pm0.5 with initial perturbations of amplitude 5p/p = 0.005 
and 6p/ p = 0.001. Note that these amplitudes lead to comparable {Sp/p = 0.005) or smaller 
{Sp/p = 0.001) amplitude modes in Fourier space compared to the 6p/p = 0.01 initialized 
modes at the higher resolution. Neither of the smaller amplitude YLRe3200Pm2 simula- 
tions showed any growth (as of 20-30 orbits), but both YLRe6400Pm0.5 simulations showed 
growth to saturation. 

From these experiments it seems that the effects of viscosity and resistivity are compa- 
rable and that the transition region between decay and growth to turbulence lies between 
Reynolds numbers of 3200 and 6400 for Pm near unity. This corresponds to a critical vertical 
wavelength, defined in terms of equation ffTTj) . of Xc/H ~ 27c /Rm^^'^ = 0.111 and 0.079, re- 
spectively. As viscosity (resistivity) is increased, MRI growth can be achieved by decreasing 
the resistivity (viscosity). This trend only works up to certain limits; if either the viscosity 
or resistivity is large enough, MRI growth is completely quenched, independent of the value 
of the other dissipation term. 



4.2. The Nonlinear Regime 

Of potentially greater interest than the linear MRI regime is the effect of viscosity and 
resistivity on fully developed MRI-driven turbulence. To study this nonlinear regime, we 
begin with model YLRe25600Pm4, a simulation with Re = 25600 and Pm = 4 at A^^ = 128, 
Ny = 200, Nz = 128 (Table [2]) that was run to 59 orbits in time. The MRI grows and the 
flow becomes fully turbulent. Averaging from t = 15 to 59 orbits gives a stress value of 
a = 0.05. We use this simulation at t = 36 orbits to initialize a series of simulations with 
different values of Re and Rm. These runs are labelled YN for y-field, nonlinear regime, and 
they are all run to 200 orbits, except for simulation YNRel2800Pm0.25, which was run to 
100 orbits. All the YN simulations are listed in Table [31 

When evolving onward from orbit 36 with modified dissipation terms, a simulation shows 
a rapid readjustment followed by either sustained turbulence at a new amplitude or decay to 
smooth flow, depending on the new values of Re and Rm. The column labeled "Turbulence?" 
in Table [3] states whether or not the given simulation has sustained turbulence. Note that 
for Rm > 1600, the turbulence is sustained except for the relatively viscous Re = 400 model. 
This critical Rm value is below the critical value obtained above for sustained growth in the 
linear regime when the resistivity and viscosity are comparable but near the critical Rm value 
in the linear regime in the absence of explicit viscosity. For simulations where turbulence is 
sustained, the column labeled "a" gives the time- and volume-averaged dimensionless stress, 
where the time average is calculated onward from orbit 50. 
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Table 2. Toroidal Field Simulations Initialized from Linear Perturbations 



Label 


Re 




Rm 


A 


= 64 


= 128 


YLRc800Pm0.5 


800 


0.5 


400 


8 


No 




YLReSOOPml 


800 


1 


800 


16 


No 




YLRe800Pm2 


800 


2 


1600 


32 


No 




YLRe800Pm4 


800 


4 


3200 


64 


No 




YLReSOOPmS 


800 


8 


6400 


128 


No 




YLRcl600Pm0.5 


1600 


0.5 


800 


16 


No 




YLRcieOOPml 


1600 


1 


1600 


32 


No 




YLRcl600Pm2 


1600 


2 


3200 


64 


No 




YLRel600Pm4 


1600 


4 


6400 


128 


No 




YLRel600Pm8 


1600 


8 


12800 


256 


No 




YLRe3200Pin0.5 


3200 


0.5 


1600 


32 


No 


No 


YLRe3200Pml 


3200 


1 


3200 


64 


No 


No 


YLRc3200Pm2 


3200 


2 


6400 


128 


Yes 


No 


YLRc3200Pm4 


3200 


4 


12800 


256 




Yes 


YLRe6400Pm0.5 


6400 


0.5 


3200 


64 


Yes 


No 


YLRe6400Pml 


6400 


1 


6400 


128 


Yes 


Yes 


YLRe6400Pm2 


6400 


2 


12800 


256 


Yes 


Yes 


YLRe6400Pm4 


6400 


4 


25600 


512 


Yes 


Yes 


YLRel2800Pm0.5 


12800 


0.5 


6400 


128 


Yes 


Yes 


YLRcl2800Pml 


12800 


1 


12800 


256 


Yes 


Yes 


YLRcl2800Pm2 


12800 


2 


25600 


512 


Yes 


Yes 


YLRcl2800Pm4 


12800 


4 


51200 


1024 


Yes 


Yes 


YLRe25600Pm0.5 


25600 


0.5 


12800 


256 


Yes 


Yes 


YLRe25600Pml 


25600 


1 


25600 


512 


Yes 


Yes 


YLRc25600Pm2 


25600 


2 


51200 


1024 


Yes 


Yes 


YLRe25600Pm4 


25600 


4 


102400 


2048 


Yes 


Yes 



The column labeled "((A))" gives a time- and volume-averaged A value in the final state 
of each simulation. Unlike Rm, A will change with the evolving magnetic field strength. 
Beginning with equation ([TU]) . we write 



(12) 



to give 



(A) 



Rm {B^) 



(13) 



where the angled brackets denote a volume average. One could also volume-average the 
square of the Alfven speed in the calculation of /? instead of averaging B'^ and p separately 
(e.g., /3 = 2cl/{v\)). We have calculated (A) using both types of averages for several frames 
in the saturated state of a few simulations. We have found at most a factor of 2 difference 
between the different calculations. Since (i?^) varies by a similar factor in the saturated 
state (see Fig. H]), this factor of 2 difference is within the uncertainty of A at any given 
time. The time- average of the volume- averaged Elsasser number, ((A)), as given in the 
table, is calculated from orbit 50 until the end of the simulation. For the decayed turbulence 
simulations in which the turbulence has not fully decayed by orbit 50, the time average is 
calculated onward from a point at which the volume-averaged magnetic energy is constant 
in time. Note that for these decayed turbulence simulations, ((A)) should equal the j3 = 100 
value associated with the net toroidal field, as given in Table |2l However, because of the 
evolution of the net By (see §[2]), the value of ((A)) after the turbulence has decayed will be 
slightly different than the j3 = 100 value. 

Since the magnetic field varies within the domain, the local value of A can also vary 
from the overall average. Histograms showing the number of grid zones with v\ of a certain 
value reveal that the percentage of grid zones that have A < 1 is at most ~ 0.01%. For the 
sustained turbulence models, ((A)) is typically on the order of 100-1000; the smallest value 
for a run with sustained turbulence is 106, and the largest value associated with a run that 
decays is 30. 

The behavior of the MRI is often characterized by the vertical component of the Alfven 
speed, and as such, we have also calculated the Elsasser number using only the vertical 
component of the magnetic field. 



(A.) 



Rm {Bfi 



(14) 
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Table 3. Toroidal Field Simulations Initialized from Nonlinear Turbulence 



Label 


Re 




Rm 


Turbulence? 


a 


((A» 




YNRc400Pm0.5 


400 


0.5 


200 


No 




4 




YNRe400Pml 


400 


1 


400 


No 




8 




YNRe400Pm2 


400 


2 


800 


No 




15 




YNRe400Pm4 


400 


4 


1600 


No 




30 




YNRc400Pm8 


400 


8 


3200 


Yes 


0.043 


614 


16.8 


YNRc400Pml6 


400 


16 


6400 


Yes 


0.068 


1983 


58.2 


YNRc800Pm0.25 


800 


0.25 


200 


No 




4 




YNRe800Pm0.5 


800 


0.5 


400 


No 




8 




YNReSOOPml 


800 


1 


800 


No 




15 




YNRe800Pm2 


800 


2 


1600 


Yes 


0.019 


137 


3.87 


YNRe800Pm4 


800 


4 


3200 


Yes 


0.038 


495 


18.0 


YNRc800Pm8 


800 


8 


6400 


Yes 


0.054 


1413 


56.2 


YNRcl600Pm0.5 


1600 


0.5 


800 


No 




15 




YNRcieOOPml 


1600 


1 


1600 


Yes 


0.018 


120 


4.45 


YNRel600Pm2 


1600 


2 


3200 


Yes 


0.033 


403 


18.6 


YNRel600Pm4 


1600 


4 


6400 


Yes 


0.044 


1120 


52.6 


YNRe3200Pm0.5 


3200 


0.5 


1600 


Yes 


0.016 


106 


4.53 


YNRe3200Pml 


3200 


1 


3200 


Yes 


0.025 


314 


16.4 


YNRc3200Pm2 


3200 


2 


6400 


Yes 


0.035 


860 


47.4 


YNRc3200Pm4 


3200 


4 


12800 


Yes 


0.043 


2170 


127 


YNRc6400Pm0.5 


6400 


0.5 


3200 


Yes 


0.021 


263 


14.9 


YNRe6400Pml 


6400 


1 


6400 


Yes 


0.031 


748 


45.2 


YNRe6400Pm2 


6400 


2 


12800 


Yes 


0.038 


1880 


118 


YNRel2800Pm0.25 


12800 


0.25 


3200 


Yes 


0.021 


262 


15.8 
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where the angled brackets denote a volume average. We have calculated the time average of 
this number, ((A^)), onward from orbit 50 for all the sustained turbulence YN simulations. 
This number is displayed in the last column of Table [31 The decayed turbulence simulations 
have approaching zero, and we do not calculate a vertical Elsasser number for these. 
Again, we calculated the vertical Elsasser number both by averaging and p separately as 
well as by averaging the ratio Bl/p. We compared the two calculations for several frames 
and found at most a factor of 1.3 difference between them. 

The ((A^)) values for the runs that have Rm closest to the critical value are on the order 



unity, with the smallest value being 3.87. As touched upon by [Fleming et al.l (120001 ) . growth 
of the vertical field MRI is largely suppressed for v\^/{rjVL) < 1 (i.e., for vertical Elsasser 
numbers less than unity). That we find ((A^)) ~ 1 close to the "decayed turbulence" regime 
may suggest that the vertical field MRI plays an important role in the sustained nonlinear 
turbulence of these toroidal field simulations. One trend to note from these data is that 
the ratio of ((A^)) to ((A)) increases with both decreasing u and decreasing 77; the vertical 
magnetic energy becomes a larger fraction of the total magnetic energy as either dissipation 
term is reduced. 

The evolution of the magnetic energy in a typical set of simulations is shown in Fig. |H 
For these runs. Re = 1600 and Rm varies by factors of two from Rm = 800 to 6400. The 
black line shows the initial evolution of YLRe25600Pm4, whose state at 36 orbits serves as 
the initial condition. It is clear that decreasing the resistivity (increasing the number) 
enhances the saturation level, and for a large enough resistivity, the turbulence decays. 

To quantify the dependence of the saturation amplitude on the dissipation coefficients, 
we plot the a values for the ensemble of simulations as a function of Re, Rm and Pm- 
Figure \5\ shows a versus Rm; the color indicates Re value, and the symbols correspond to 
the Pm value. The simulations with a = are those where the turbulence decayed away, 
which include all simulations with Rm < 800 and the Re = 400, Rm = 1600 simulation. 
Overall there is a general trend of increasing a value with decreasing resistivity. 

The dependence of a on Re is shown in Fig. O Here the color indicates the Rm value, 
whereas Pm is again represented by a symbol. Evidently, if the resistivity is low enough, 
increasing the viscosity will increase the a values. However, consider the YN simulations 
with Rm = 1600. These simulations suggest that if the resistivity is close to some critical 
value, increasing the viscosity will cause the turbulence to decay. Another feature of note is 
that as Re increases, the range of a for different Rm values becomes smaller, and a appears 
to converge to ~ 0.02 — 0.04 for all Rm. This could indicate that as u and 77 decrease, their 
influence on the turbulence level might decrease. However, for large values of Re or Rm, 
the dissipation lengthscales are under-resolved, and higher resolution is needed to test this 
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possibility. 



We plot the dependence of a on P^n in Fig. [71 In this figure, the colors represent varying 
Rm while the symbols denote different Re values. The clearest trend is that if Rm is large 
enough to sustain turbulence, increasing Pm leads to larger a values. Note that turbulence 
can be sustained even for Pm less than unity, if Rm is large enough. At constant Pm, we 
find that a oc Pe-^i with 5i ranging from -0.1 to -0.3 (calculated by a linear fit to the data 
in log space for non-decayed turbulence simulations only). At constant Re value, we find 
OL oc Rm^"^ with in the range 0.4-0.8 and generally decreasing with increasing Re. 

These results naturally lead to the question of why increasing v or decreasing t] causes 
an increase in turbulence. Magnetic reconnection and dissipation of field lines, either due to 
an explicit resistivity or to grid-scale effects, presumably play the primary role in limiting 
the amplitude of the MHD turbulence. iBalbus fc Hawleyl (Il998l ) hypothesized that increased 
viscosity would inhibit reconnection by preventing velocity motions that would bring field 
together on small scales. When P^ > 1, the viscous length is greater than the resistive 
one, and m agnetic field dissipatio n becomes less efficient, leading to an increase in turbulent 
stress (e.g.. iBalbus fc Henrill2008l ). If this hypothesis is correct, there may also be a change 
in the dissipation of kinetic and magnetic energy into heat. To investigate this possibility, 
we carry out an analysis of viscous and resistive heating for several of the simulations. 



Consider the vol ume-averaged 



tions (15) and (16) in lSimon et al 



dneti c and magnetic energy evolution equations, equa- 
J2009h . 



K 



- V- 



P + 1 V ■ V 



B[v B 

B-{B-Vv))-G-Q^ 



(15) 



and 



M = -( V 



-B'v 



B^V-v) + {B-{B-Vv))-Q, 



(16) 



Here, K and M are the time derivatives of the volume-averaged kinetic and magnetic ener- 
gies, respectively. The time derivative of the volume-averaged gravitational potential energy 
is given by G, and Qk and Qm are the volume-averaged kinetic and magnetic energy dissi- 
pation rates, respectively. The gravitational potential is $ = qQ^{^ 



X 
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We determine Q\ and Qm for select YN models by computing the time average of each 
of the source terms in the energy evolution equations using 200 data files equa lly spaced in 



time over 20 orbits. We assume G is zero in the time-average; the analysis of ISimon et al. 



(120091 ) found G is always negligibly small. The time-derivatives, K and M, are calculated by 
differentiating the volume-averaged kinetic and magnetic energy history data with respect 
to time and then sampling these data to the times associated with the data files. The 
dissipation terms Qk and Qm, which include both physical and numerical effects, are the 
remainder after all the other terms are calculated. 

Figure [8] shows the ratio of the time-average (Qk) to (Qm) as a function of Pm and a 
for select YN runs. The colors and symbols are the same as in Fig. [61 The time average is 
calculated from t = 70-90 orbits for YNRe400Pml6 (black X) and YNRel2800Pm0.25 (blue 
circle), t = 110 - 130 orbits for YNRe800Pm2 (green diamond) and YNReSOOPmS (black 
square), and t = 110.6 - 130.6 orbits for YNRe800Pm4 (blue triangle) and YNRe3200Pm4 
(red triangle). The ratio of viscous to resistive heating generally increases as either a or Pm 
increases, although not monotonically. The relative heating ratio is not simply proportional 
to Pm as one might naively expect. 



The data suggest a general relationship between saturated stress a nd (Qy) / (Q 



We 



know that the stress level sets the total dissipation rate (Qk + Qm) (e.g.. lSimon et al.ll2009l ): 
stronger stresses extract more energy from the background shear flow, and that turbulence is 
rapidly dissipated into heat. However, does stronger turbulence by itself change the heating 
ratio, or is the change in the heating ratio mainly determined by changes in Pm, which also 
increase the turbulence levels? This question of causality cannot be definitively answered 
from these data. 

Further insight may come from examining the ratio of averaged Reynolds stress, {{pVxSvy)), 
to averaged Maxwell stress, {{—B^By)), as a function of a; this is shown in Fig. O The colors 
and symbols are the same as in Fig. [HI The double brackets for the stresses denote time and 
volume averages, where the time average is calculated over the same 20 orbit period as in 
Fig. [HI There is a decrease in the ratio of the Reynolds to Maxwell stress as the total stress 
increases. These stresses are proportional to the perturbed magnetic and kinetic energies 
at the largest scales, and if this continued down to the dissipation scale, we might expect 
that the ratio (Qk) / (Qm) would behave similarly with a. In fact, the heating ratio shows 
the opposite trend with a, indicating that a transfer of energy from magnetic to kinetic 
fluctuations must occur in the turbulent cascade. 

Past net toroidal field simulations without explicit dissipation terms also fi nd a trend for 



a dec rease in the ratio of the Reynolds to Maxwell stress with increasing a (e.g.. lHawley et al. 



19951 ). So this may be a general result independent of Pm. The quantity (Qk) / {Qm) has not 
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been extensively studied in past shearing box simulations, but ISimon et al.l (120091 ) found a 



ratio of ~ 0.6 for a net vertical field model without explicit dissipation terms. 

In summary, these observations are consistent with the hypothesis that decreasing 
increases the efficiency of magnetic reconnection and hence reduces the overall stress level. 
However, a more in-depth study would be required to better understand the full causal 
relationship between the ratio of dissipation terms and the saturation levels. 

Finally, we note that the ratio of Reynolds stress to perturbed kinetic energy increases 
with increasing z/, as shown in Fig. [TOl There is no observed trend with rj. As u is increased, 
the fluid motions that are not being directly driven by the MRI become increasingly damped. 
The fluid motions that are driven by the magnetic field in the form of Reynolds stress follow 
the behavior of the Maxwell stress with u. This is also consistent with the hypothesis that 
increased u leads to less efficient magnetic reconnection; the kinetic fluctuations become 
damped relative to the driving via the MRI, making it difficult to bring field lines close 
together for reconnection. 

Overall, resistivity seems to play a more fundamental role than viscosity in these net 
toroidal field simulations. There is a critical Rm below which turbulence decays or fails to 
grow from linear perturbations. For a given resistivity near this critical value, a relatively low 
viscosity leads to MRI growth (linear regime) or sustained turbulence (nonlinear regime). A 
high viscosity can prevent growth (linear regime) or cause decay (nonlinear regime). Once the 
resistivity is sufficiently low to ensure MRI growth to saturation and continued turbulence, 
the effect of viscosity changes and higher viscosity gives larger a values. 



5. Discussion and Conclusions 

In this study, we carried out a series of local, unstratified shearing box simulations of the 
MRI with Athena including the effects of constant shear viscosity and Ohmic resistivity. The 
first simulations were initialized with a zero net magnetic fiux in the domain for comparison 
with the results of F07. The second set of simulations are the first investigation of the impact 
of both viscosity and resistivity on models with a net toroidal field. 

For the values of viscosity and resistivity they studied, F07 found that turbulence was 
sustained only above a critical value, specifically when > 1. There was evidence 
that this critical Pm value decreases as Re increases (viscosity is reduced). We repeated 
these experiments and found that the saturation level of the MRI depends strongly on both 
viscosity and resistivity, and for every Re, there exists a critical Pm value below which the 
turbulence dies out. For those simulations where turbulence was sustained, we found good 
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agreement between the Athena a values and those of F07. 

Zero net field simulations are fundamentally different from net field models because 
an imposed background field cannot be removed as a result of the simulation evolution. 
The net field remains unstable and can drive fiuid motion even during the fully nonlinear 



turbulence phase, assuming that that field was unstable to begin with. iLesur &: Longaretti 



(120071 ) examined the effects of diffusion on models with a net /3 = 100 vertical field in a 
1x4x1 shearing box using a pseudo-spectral incompressible code. They found a relation 
a cx with b = 0.25-0.5 for values of ranging from 0.12 to 8, but they found no case 
where the turbulence died out completely for the range of viscosities and resistivities studied. 

Among net field models, the purely vertical field case is significantly different from the 
purely toroidal field model, hence the need for the study we have presented here. For a 
vertical field, the linear MRI favors wavenumbers ~ ^/va and = ky = 0. The purely 
toroidal case favors ky ~ Q/va with k/kz minimized. Since k^ is time dependent due to the 
background shear, a given mode undergoes a finite period of amplification as kx swings from 
leading to traihng. These properties suggest that purely toroidal field models might be more 
sensitive to dissipation than the vertical field case. 

In our numerical study of the linear growth regime of the toroidal MRI, we have found 
that increasing either the viscosity or the resistivity can prevent the growth of MRI modes. 
As the viscosity (resistivity) increases, the MRI needs a smaller resistivity (viscosity) in 
order to grow. However, for large enough values of either the viscosity or the resistivity, 
MRI growth is suppressed, even in the absence of the other dissipation term. Because 
of the importance of small wavelength (large wavenumber) modes, the critical Rm values, 
below which growth is inhibited, tends to be larger than what one would expect from an 
axisymmetric vertical field analysis, even in the absence of viscosity. Here, for comparable 
values of viscosity and resistivity, the critical Rm value was around 3200-6400. 

Because the linear toroidal field MRI is time dependent, turbulence can only be sustained 
if nonlinear amplitudes are reached during the growth phase. Thus, the outcome of the linear 
MRI phase can be sensitive to the initial amplitude of the perturbations in a simulation where 
the viscous or resistive values are near the critical value. 

In the nonlinear regime, we found that viscosity generally acts in an opposite sense 
to that in the linear regime; increased viscosity enhances angular momentum transport. 
Furthermore, increasing the resistivity appears to decrease the saturation level, in agreement 
with previous studies, and the critical Rm, below which the turbulence dies is ~ 800-1600. 
Near the critical Rm, however, increasing the viscosity causes the turbulence to decay, a 
behavior more in line with the linear regime. 
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In our simulations, as well as those of iLesur fc Longarettil (120071 ). Pm < 1 did not 



necessarily quench the nonlinear turbulence or prevent growth from linear perturbations. 
Resistivity or viscosity above a certain level can stabilize the system against these pertur- 
bations, but if both are sufficiently small, their ratio has no influence on MRI growth. The 
presence of turbulence, however, is distinct from the saturation level of that turbulence, and 
here Pm can have a significant effect. For both net toroidal and net vertical field simulations, 
a increases with increasing P^ for the range of values studied. 

What do these results imply for the effect of resistivity and viscosity on the MRI and 
on astrophysical systems? In principle, they could be quite significant. In protostellar disks, 
the resistiv ity is thought to be quite high near the midplane, leading to the existence of the 



dead zone (lGammielll996l ). The Rm values studied here could be applicable to such systems. 
However, the implications for accretion disks with small values of viscosity and resistivity 
(e.g.. X-ray binary disks) are less clear. Because the range of a values we obtained decreases 
with increasing Re (Fig. [6]), it is possible that a may converge to a single value independent 
of Pm as Re and Rm are increased. If true, this would suggest that the dissipation scales 
might have very little infiuence on the saturation level of the MRI in astrophysical disks. 
This idea will need to be tested with higher resolution simulations to ensure that the (small) 
viscous and resistive scales are adequately resolved. If, on the other hand, Pm still has an 
infiuence on the turbulence even for very small values of viscosity and resistivity, our results 
(taken together with those in the literature) could be applicable to such disks. The resistivity, 
viscosity, and Pm can vary quite substantially in these systems, not only between differen t 



astrophysical ob j ects, b ut also within a given disk (e.g.. lBrandenburg &: Subramanianll2005l ). 



Balbus fc Henri! (120081 ) analyze a possible Pm dependence on radius in X-ray binaries to 
suggest that such a dependence could be at the core of spectral state transitions in these 
systems. 

The work to date is suggestive, but there remain several limitations associated with 
these shearing box simulations. First, the simulations are unstratified; v ertical gravi t y may 



also play a role in establishing the overall turbulent state. For example, iDavis et al.l (120091 ) 
carried out a series of zero net field shearing box simulations with vertical gravity and explicit 
dissipation and found that the turbulence does not decay as readily as in the unstratified 
case. Secondly, all of the simulations to date have explored a relatively restricted range 
of parameters. Here, for example, we have examined only one value for the toroidal field 
strength and one domain size. Finally, as touched upon above, the range of values for Re 
and Rm that have been studied are somewhat restricted and certainly much smaller than 
would be appropriate for many astrophysical disk systems. While this limitation is partially 
computational and can be improved upon with higher resolutions, the question remains for 
astrophysical systems whether viscous and resistive processes that take place on relatively 
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small lengthscales can have a significant influence on macroscopic stress terms whose scales 
are on order the pressure scale height in the disk. But regardless of the importance of resis- 
tivity and viscosity for astrophysical systems, the values of Re and Rm are very important 
for setting a in numerical simulations, much more so than many other shearing box param- 
eters (e.g., pressure) studied to date. Without a more thorough understanding of the role 
that dissipation terms play, quantitative predictions of a values from simulations will not be 
possible. 

In summary, our experiments have explored the effect of changing viscosity and resistiv- 
ity on MRl simulations with a net toroidal field. This work serves to expand upon previous 
investigations of the impact of small-scale dissipation. While the direct applicability of stud- 
ies such as this to specific stress values within astrophysical systems remains uncertain, it is 
likely that for the conceivable future, numerical simulations will be our primary, if not only 
way to explore the nature of MRI-driven turbulence. A thorough understanding of MRI 
turbulence can only be obtained with a complete understanding of the effects of diffusion, 
both numerical and physical. 
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Fig. 1. — Numerical error as a function of x resolution for the three-dimensional decaying 
linear wave problem. The boxes are the errors for a decaying Alfven wave, and the asterisks 
are the errors for a decaying sound wave. The error is calculated from the square root of the 
sum of the squared errors in the density and momenta (for the sound wave) and the density, 
momenta, and magnetic field (for the Alfven wave) obtained using the analytic solution. 
The dashed line shows the slope corresponding to second-order convergence. 
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Fig. 2. — Time- and volume-averaged stress parameter a as a function of Pm in the SZ 
simulations; a = {{pVxSvy — B^By)) / Po, where the average is calculated over the entire 
simulation domain and from 20 orbits to the end of the simulation. Only simulations with 
sustained turbulence are plotted. The Pm = 4 model has Rm = 12500 whereas the other 
two have Rm — 12800. There is a nearly linear increase in a with Pm. 
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Fig. 3. — Time evolution of volume-averaged magnetic energy density normalized by the gas 
pressure for the SZ runs with Re — 1600 and varying Pm- The volume average is calculated 
over the entire simulation domain. The solid line corresponds to Pm = 8, the dashed line 
corresponds to = 4, and the dotted line corresponds to Pm = 2. The turbulence decays 
for the lowest two Pm values, with the Pm — 4 case taking roughly 60 orbits to decay. 
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Fig. 4. — Time evolution of volume-averaged magnetic energy density normalized by the gas 
pressure for the YN runs with Re = 25600 (black curve) and Re = 1600 (colored curves). 
The volume average is calculated over the entire simulation domain. The colors indicate 
Pm; green corresponds to Rm = 800 (Pm = 0.5), light blue to Rm = 1600 (Pm = 1), red to 
Rm = 3200 (Pm = 2), and dark blue to Rm = 6400 (Pm = 4). Increasing Rm (Pm) leads to 
enhanced turbulence. 
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Fig. 5. — Time- and volume-averaged stress parameter a as a function of Rm in the YN 
simulations; a = {{pv^^Vy — B^By)) / Pg. The time average runs from 50 orbits onward, and 
the volume average is calculated over the entire simulation domain. The colors correspond 
to Re values, and the symbols correspond to Pm values. Red symbols are Re = 400, green 
Re = 800, dark blue Re = 1600, black Re = 3200, pink Re = 6400, and light blue are 
Re = 12800. Circles are = 0.25, crosses Pm = 0.5, asterisks Pm = 1, diamonds Pm = 2, 
triangles Pm = 4, squares Pm = 8, and X's are Pm = 16. Note that some of the decayed 
turbulence {a = 0) simulations are not plotted for clarity. Increasing Rm results in larger a 
values, and for Rm less than 800-1600, the turbulence decays. 
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Fig. 6. — Time- and volume-averaged stress parameter a as a function of Re in the YN 
simulations; a = {{pv^Svy — B^By)) / Pq. The time average runs from 50 orbits onward, and 
the volume average is calculated over the entire simulation domain. The colors correspond 
to Rm values, and the symbols correspond to Pm values. Light blue symbols are Rm = 800, 
green Rm = 1600, dark blue Rm = 3200, black Rm = 6400, and red are Rm = 12800. 
Circles are Pm = 0.25, crosses Pm = 0.5, asterisks Pm = 1, diamonds Pm = 2, triangles 
Pm = 4, squares Pm = 8, and X's are Pm = 16. Note that some of the decayed turbulence 
{a = 0) simulations are not plotted for clarity. Increasing Re leads to decreasing a values. 
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Fig. 7. — Time- and volume-averaged stress parameter a as a function of P^; a = 
{{pVxSvy — BxBy)) / Pq. The time average runs from 50 orbits onward, and the volume average 
is calculated over the entire simulation domain. The colors correspond to Rm values, and 
the symbols to Re values. Light blue symbols are Rm = 800, green Rm = 1600, dark blue 
Rm = 3200, black Rm = 6400, and red are Rm = 12800. Crosses are Re = 400, asterisks 
Re = 800, diamonds Re = 1600, triangles Re = 3200, squares Re = 6400, and circles are 
Re = 12800. Note that some of the decayed turbulence {a = 0) simulations are not plotted 
for clarity. The average stress increases with increasing P^. 
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Fig. 8. — Ratio of kinetic to magnetic energy dissipation rate as a function of Pm (top panel) 
and a (bottom panel) for select YN runs; a = {{pv^Svy — B^By)) / Pq. The colors and symbols 
are the same as in Fig. [61 The kinetic and magnetic dissipation rates as well as a have been 
averaged in volume and time. The volume average is calculated over the entire simulation 
domain and the time average is calculated from t = 70 — 90 orbits for YNRe400Pml6 (black 
X) and YNRel2800Pm0.25 (blue circle), t = 110-130 orbits for YNRe800Pm2 (green 
diamond) and YNRe800Pm8 (black square), and t = 110.6 - 130.6 orbits for YNRe800Pm4 
(blue triangle) and YNRe3200Pm4 (red triangle). The ratio of viscous to resistive heating 
generally increases as either a or increases, although not monotonically. 
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Fig. 9. — Ratio of Reynolds stress to Maxwell stress as a function of a for select YN runs; 
a = {{pVxSvy — B^By)) / Po- The colors and symbols are the same as in Fig. O The Maxwell 
and Reynolds stresses as well as a have been averaged in volume and time. The volume 
average is calculated over the entire simulation domain and the time average is calculated 
from t = 70 - 90 orbits for YNRe400Pml6 (black X) and YNRel2800Pm0.25 (blue circle), 
t = 110 - 130 orbits for YNRe800Pm2 (green diamond) and YNReSOOPmS (black square), 
and t = 110.6 - 130.6 orbits for YNRe800Pm4 (blue triangle) and YNRe3200Pm4 (red 
triangle). The ratio of Reynolds to Maxwell stress generally decreases with increasing a. 
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Fig. 10. — Ratio of Reynolds stress to perturbed kinetic energy as a function of Re in the 
sustained turbulence YN simulations. Both the Reynolds stress and the perturbed kinetic 
energy are time and volume averaged, with the time average calculated from orbit 50 onward 
and the volume average calculated over the entire simulation domain. The colors correspond 
to Rm values, and the symbols correspond to Pm values. Green symbols are Rm = 1600, 
blue Rm = 3200, black Rm = 6400, and red are Rm = 12800. Circles are Pm = 0.25, crosses 
Pm = 0.5, asterisks Pm = 1, diamonds Pm = 2, triangles P^ = 4, squares Pm = 8, and X's 
are Pm = 16. The ratio of stress to energy increases with increasing viscosity but does not 
change with resistivity. 



